
regressIter = function(outcome, iv1, iv2=NA, controls, hawkType="hawkBoot", theFamily, dataset) {
  
  maindata = dataset |> filter(hawkMeasure==all_of(iv1))
  
  if (is.na(iv2)) {
    
    # Iterate
    alliters = foreach (i = 1:1000, .combine="rbind") %dopar% {
      
      # Get the one hawkishness column
      onecol = paste0(hawkType, i)
      
      # Get the data with the hawkishness and DV columns properly named
      subdata = maindata |> dplyr::select(dv=all_of(outcome), hawk1=all_of(onecol), everything()) |> 
        dplyr::select(-starts_with("hawkBoot")) 
      names(subdata)[which(names(subdata)=="hawk1")] = iv1
      
      # Get the formula
      theFormula = paste0("dv", "~", iv1, paste(controls, sep="+"))
      
      # Do the regression
      onefit = glm(theFormula, data=subdata, family=theFamily)
      
      # Store all results in a list
      results = broom::tidy(onefit)
      results$nobs = nobs(onefit)
      
      results
    }
  } else {
    data2 = dataset |> filter(hawkMeasure==all_of(iv2)) |> dplyr::select(starts_with("hawkBoot"))
    names(data2) = paste0(names(data2), "iv2")
    
    maindata = bind_cols(maindata, data2)
    
    # Iterate
    alliters = foreach (i = 1:1000, .combine="rbind") %dopar% {
      
      # Get the one hawkishness column
      onecol = paste0("hawkBoot", i)
      twocol = paste0("hawkBoot", i, "iv2")
      
      # Get the data with the hawkishness and DV columns properly named
      subdata = maindata |> dplyr::select(dv=all_of(outcome), hawk1=all_of(onecol), hawk2=all_of(twocol), everything()) |> 
        dplyr::select(-starts_with("hawkBoot"))
      names(subdata)[which(names(subdata)=="hawk1")] = iv1
      names(subdata)[which(names(subdata)=="hawk2")] = iv2
      
      # Get the formula
      theFormula = paste0("dv", "~", iv1, "+", iv2, paste(controls, sep="+"))
      
      # Do the regression
      onefit = glm(theFormula, data=subdata, family=theFamily)
      
      # Store all results in a list
      results = broom::tidy(onefit)
      results$nobs = nobs(onefit)
      results
    }
  }

  # Record number of observations in the analysis (same across all models of same class)
  nobservations = alliters$nobs[1]
  
  # Put together all results
  alliters_sub = alliters |> filter(!grepl("mention|admin", term))
  
  # Relabel and reorganize variables of interest
  alliters_sub$term = factor(alliters_sub$term,
                             levels=c("meanHawk", "advHawk", "presHawk", 
                                      "meanHawkKey", "advHawkKey", "presHawkKey", 
                                      "nAttend", 
                                      "numDef", "numInt", "numMil", "numState",
                                      "log(totDip + 1)", "log(totInt + 1)", "log(totMil + 1)",
                                      "log(usmidschallenge5 + 1)", "cinc", "formal",
                                      "factor(admin)Eisenhower", "factor(admin)Kennedy", "factor(admin)Johnson", "factor(admin)Nixon",
                                      "factor(admin)Ford", "factor(admin)Carter", "factor(admin)Reagan", "(Intercept)"),
                             labels=c("Mean Hawkishness", "Advisers' Hawkishness (Acts)", "President's Hawkishness", 
                                      "Mean Hawkishness", "Advisers' Hawkishness (Acts)", "President's Hawkishness", 
                                      "No. of Attendees", 
                                      "Defense", "Intelligence", "Military", "State",
                                      "Diplomatic Experience", "Intelligence Experience", "Military Experience",
                                      "5-Year MID Challenges", "US CINC", "Formal", 
                                      "Eisenhower", "Kennedy", "Johnson", "Nixon",
                                      "Ford", "Carter", "Reagan",
                                      "Constant"))
  
  # Get key terms to summarize
  terms = unique(alliters_sub$term)
  
  # Summarize
  termList = list()
  for (j in 1:length(terms)) {
    oneterm = alliters_sub |> filter(term==terms[j]) |> dplyr::select(estimate, std.error)
    
    termstats = oneterm |> summarize(est_mean = mean(estimate, na.rm=T), est_lower = quantile(estimate, 0.025, na.rm=T), est_upper = quantile(estimate, 0.975, na.rm=T),
                                     se_mean = mean(std.error, na.rm=T), se_lower = quantile(std.error, 0.025, na.rm=T), se_upper = quantile(std.error, 0.975, na.rm=T))
    
    termList[[j]] = data.frame(term=terms[j], termstats)  
  }
  allTerms = rbindlist(termList) 
  
  
  allOutput = list(results=allTerms, nobs=nobservations)

  return(allOutput)
}




coefTable = function(onemodelresult) {
  
  # Get the terms in the regression, in desired order
  modelTerms = onemodelresult[["results"]][["term"]]
  
  terms = factor(modelTerms,
                 levels=c("Mean Hawkishness", "Advisers' Hawkishness (Acts)", "President's Hawkishness", 
                          "No. of Attendees", 
                          "Defense", "Intelligence", "Military", "State",
                          "Diplomatic Experience", "Intelligence Experience", "Military Experience",
                          "5-Year MID Challenges", "US CINC", "Formal", 
                          "Eisenhower", "Kennedy", "Johnson", "Nixon",
                          "Ford", "Carter", "Reagan",
                          "Constant"))
  terms = levels(terms)[which(levels(terms) %in% modelTerms)]
  
  # Get the core statistics and significance
  ceList = list()
  for (i in 1:length(terms)) {
    oneterm = onemodelresult[["results"]] |> filter(term==terms[i])
    
    onetermest = sprintf("%.3f", round(oneterm$est_mean, 3))
    onetermse = sprintf("%.3f", round(oneterm$se_mean, 3))
    
    # Add another digit if numbers are too small
    onetermest = ifelse(onetermest=="0.000", sprintf("%.4f", round(oneterm$est_mean, 4)), onetermest)
    onetermse = ifelse(onetermse=="0.000", sprintf("%.4f", round(oneterm$se_mean, 4)), onetermse)
    
    # Put parentheses around the standard errors
    onetermse = paste0("(", onetermse, ")")
    
    # Find significance 
    onelower90 = oneterm$est_mean - qnorm(0.95)*oneterm$se_mean
    oneupper90 = oneterm$est_mean + qnorm(0.95)*oneterm$se_mean  
    
    onelower95 = oneterm$est_mean - qnorm(0.975)*oneterm$se_mean
    oneupper95 = oneterm$est_mean + qnorm(0.975)*oneterm$se_mean  
    
    onelower99 = oneterm$est_mean - qnorm(0.995)*oneterm$se_mean
    oneupper99 = oneterm$est_mean + qnorm(0.995)*oneterm$se_mean  
    
    stars = ifelse(onelower99 > 0 | oneupper99 < 0, "^{***}",
                   ifelse(onelower95 > 0 | oneupper95 < 0, "^{**}",
                          ifelse(onelower90 > 0 | oneupper90 < 0, "^{*}", "")))
    onetermest_star = paste0(onetermest, stars) 
    
    # Put things together for one term  
    onetermres = data.frame(rbind(t(onetermest_star), t(onetermse)))
    onetermres$varToMatch = c(paste0(terms[i], "_est"), paste0(terms[i], "_se"))
    onetermres = onetermres |> dplyr::select(varToMatch, everything())
    names(onetermres) = c("varToMatch", "est") 
    
    ceList[[i]] = onetermres
    
  }
  ceOut = rbindlist(ceList)
  
  ceOut = rbind(ceOut, data.frame(varToMatch="Observations", est=onemodelresult$nobs))
  
  # Put out the results table
  return(ceOut)
  
  # Print key results
  #print(xtable(ceOut), sanitize.text.function=identity, include.rownames=F) 
  
}


makeFullCoefTable = function(coefList) {
  
  # How many tables need to be combined?
  nTables = length(coefList)

  # Get first coefficient table
  basetable = coefTable(coefList[[1]])
  
  # Attach additional tables to the base line
  for (i in 2:nTables) {
    anothertable = coefTable(coefList[[i]])
    
    basetable = merge(basetable, anothertable, by="varToMatch", all=T)
    names(basetable) = c("varToMatch", paste0("est", 1:i))
  }

    
  # Need to reorder variables
  basetable$var = gsub("_est|_se", "", basetable$varToMatch)
  basetable$id = NA
  basetable$id[grepl("_est", basetable$varToMatch)] = 1
  basetable$id[grepl("_se", basetable$varToMatch)] = 2
  basetable$id[grepl("Observations", basetable$varToMatch)] = 3
  
  basetable$var = factor(basetable$var,
                 levels=c("Mean Hawkishness", "Advisers' Hawkishness (Acts)", "President's Hawkishness", 
                          "No. of Attendees", 
                          "Defense", "Intelligence", "Military", "State",
                          "Diplomatic Experience", "Intelligence Experience", "Military Experience",
                          "5-Year MID Challenges", "US CINC", "Formal", 
                          "Eisenhower", "Kennedy", "Johnson", "Nixon",
                          "Ford", "Carter", "Reagan",
                          "Constant", "Observations"))
  
  basetable = basetable |> arrange(var, id) |> dplyr::select(var, varToMatch, starts_with("est")) 
  basetable$var[grepl("_se", basetable$varToMatch)] = ""
  basetable = basetable |> dplyr::select(-varToMatch)
  
  print(xtable(basetable), include.rownames=F, type='latex', sanitize.text.function=identity)
  #return(basetable)
}



######### MONTHLY ANALYSIS #########
monthIter = function(controls, hawkType="hawkBoot", monthlydata=monthdata, presdata = presHawkData, advdata = hmcData) {
  
  # Iterate
  alliters = foreach (i = 1:1000, .combine="rbind") %dopar% {
    
    # Get the one hawkishness column
    onecol = paste0(hawkType, i)
    
    presColumn = presdata |> dplyr::select(all_of(onecol)) |> unlist()
    advColumn = advdata |> dplyr::select(all_of(onecol)) |> unlist()
    
    monthdata2 = data.frame(monthlydata, presHawk=presColumn, mainHawk=advColumn)
    
    theFormula = paste0("newmids_us ~ mainHawk + presHawk", controls)
    
    midDep = glm(theFormula, data=monthdata2, family=poisson)
    
    # Store all results in a list
    results = broom::tidy(midDep)
    results$nobs = nobs(midDep)
    results
  }
  
  # Record number of observations in the analysis (same across all models of same class)
  nobservations = alliters$nobs[1]
  
  # Put together all results
  alliters_sub = alliters |> filter(!grepl("mention|admin", term))
  
  # Relabel and reorganize variables of interest
  alliters_sub$term = factor(alliters_sub$term,
                             levels=c("mainHawk", "presHawk",
                                      "war.now", "log.last.war.dead.pc", "log.months.since.last.war",
                                      "last.war.win", "n.mid.us.challenge.past5", "prop.mid.win.us.past5",
                                      "is.recession", "unigov", "cinc", "log(tenuremonth)", "(Intercept)"),
                             labels=c("Advisers' hawkishness", "President's hawkishness",
                                      "War ongoing", "Deaths per capita in last war (logged)", "Months since last war (logged)",
                                      "Victory in last war", "MID challenges to US in last 5 years", "Average MID outcome in last 5 years",
                                      "Economic recession", "Unified government", 
                                      "US material capabilities", "President's tenure (logged months)", "Constant"))
  
  # Get key terms to summarize
  terms = unique(alliters_sub$term)
  terms = terms[order(terms)]
  
  # Summarize
  termList = list()
  for (j in 1:length(terms)) {
    oneterm = alliters_sub |> filter(term==terms[j]) |> dplyr::select(estimate, std.error)
    
    termstats = oneterm |> summarize(est_mean = mean(estimate, na.rm=T), est_lower = quantile(estimate, 0.025, na.rm=T), est_upper = quantile(estimate, 0.975, na.rm=T),
                                     se_mean = mean(std.error, na.rm=T), se_lower = quantile(std.error, 0.025, na.rm=T), se_upper = quantile(std.error, 0.975, na.rm=T))
    
    termList[[j]] = data.frame(term=terms[j], termstats)  
  }
  allTerms = rbindlist(termList) 
  
  
  allOutput = list(results=allTerms, nobs=nobservations)
  
  return(allOutput)
}



# Single coefficient table code
monthCoefTable = function(onemodelresult) {
  
  # Get the terms in the regression, in desired order
  modelTerms = onemodelresult[["results"]][["term"]]
  
  
  terms = factor(modelTerms,
                 levels=c("hmcMain", "presHawk",
                          "war.now", "log.last.war.dead.pc", "log.months.since.last.war",
                          "last.war.win", "n.mid.us.challenge.past5", "prop.mid.win.us.past5",
                          "is.recession", "unigov", "cinc", "log(tenuremonth)", "(Intercept)"),
                 labels=c("Advisers' hawkishness", "President's hawkishness",
                          "War ongoing", "Deaths per capita in last war (logged)", "Months since last war (logged)",
                          "Victory in last war", "MID challenges to US in last 5 years", "Average MID outcome in last 5 years",
                          "Economic recession", "Unified government", 
                          "US material capabilities", "President's tenure (logged months)", "Constant"))
  
  terms = levels(terms)[which(levels(terms) %in% modelTerms)]
  
  # Get the core statistics and significance
  ceList = list()
  for (i in 1:length(terms)) {
    oneterm = onemodelresult[["results"]] |> filter(term==terms[i])
    
    onetermest = sprintf("%.3f", round(oneterm$est_mean, 3))
    onetermse = sprintf("%.3f", round(oneterm$se_mean, 3))
    
    # Add another digit if numbers are too small
    onetermest = ifelse(onetermest=="0.000", sprintf("%.4f", round(oneterm$est_mean, 4)), onetermest)
    onetermse = ifelse(onetermse=="0.000", sprintf("%.4f", round(oneterm$se_mean, 4)), onetermse)
    
    # Put parentheses around the standard errors
    onetermse = paste0("(", onetermse, ")")
    
    # Find significance 
    onelower90 = oneterm$est_mean - qnorm(0.95)*oneterm$se_mean
    oneupper90 = oneterm$est_mean + qnorm(0.95)*oneterm$se_mean  
    
    onelower95 = oneterm$est_mean - qnorm(0.975)*oneterm$se_mean
    oneupper95 = oneterm$est_mean + qnorm(0.975)*oneterm$se_mean  
    
    onelower99 = oneterm$est_mean - qnorm(0.995)*oneterm$se_mean
    oneupper99 = oneterm$est_mean + qnorm(0.995)*oneterm$se_mean  
    
    stars = ifelse(onelower99 > 0 | oneupper99 < 0, "^{***}",
                   ifelse(onelower95 > 0 | oneupper95 < 0, "^{**}",
                          ifelse(onelower90 > 0 | oneupper90 < 0, "^{*}", "")))
    onetermest_star = paste0(onetermest, stars) 
    
    # Put things together for one term  
    onetermres = data.frame(rbind(t(onetermest_star), t(onetermse)))
    onetermres$varToMatch = c(paste0(terms[i], "_est"), paste0(terms[i], "_se"))
    onetermres = onetermres |> dplyr::select(varToMatch, everything())
    names(onetermres) = c("varToMatch", "est") 
    
    ceList[[i]] = onetermres
    
  }
  ceOut = rbindlist(ceList)
  
  ceOut = rbind(ceOut, data.frame(varToMatch="Observations", est=onemodelresult$nobs))
  
  # Put out the results table
  return(ceOut)
}


# Put all results from monthly analysis together
makeMonthFullCoefTable = function(coefList) {
  
  # How many tables need to be combined?
  nTables = length(coefList)
  
  # Get first coefficient table
  basetable = monthCoefTable(coefList[[1]])
  
  # Attach additional tables to the base line
  for (i in 2:nTables) {
    anothertable = monthCoefTable(coefList[[i]])
    
    basetable = merge(basetable, anothertable, by="varToMatch", all=T)
    names(basetable) = c("varToMatch", paste0("est", 1:i))
  }
  
  
  # Need to reorder variables
  basetable$var = gsub("_est|_se", "", basetable$varToMatch)
  basetable$id = NA
  basetable$id[grepl("_est", basetable$varToMatch)] = 1
  basetable$id[grepl("_se", basetable$varToMatch)] = 2
  basetable$id[grepl("Observations", basetable$varToMatch)] = 3
  
  basetable$var = factor(basetable$var,
                         levels=c("Advisers' hawkishness", "President's hawkishness",
                                  "War ongoing", "Deaths per capita in last war (logged)", "Months since last war (logged)",
                                  "Victory in last war", "MID challenges to US in last 5 years", "Average MID outcome in last 5 years",
                                  "Economic recession", "Unified government", 
                                  "US material capabilities", "President's tenure (logged months)", "Constant", "Observations"))
  
  basetable = basetable |> arrange(var, id) |> dplyr::select(var, varToMatch, starts_with("est")) 
  basetable$var[grepl("_se", basetable$varToMatch)] = ""
  basetable = basetable |> dplyr::select(-varToMatch)
  
  print(xtable(basetable), include.rownames=F, type='latex', sanitize.text.function=identity)
  #return(basetable)
}


      
    